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ABSTRACT 

In the standard paradigm for cosmological structure formation, clustering develops 
from initially random-phase (Gaussian) density fluctuations in the early Universe by 
a process of gravitational instability. The later, non-linear stages of this process in- 
volve Fourier mode-mode interactions that result in a complex pattern of non-random 
phases. We present a novel mapping technique that reveals mode coupling induced 
by this form of nonlinear interaction and allows it to be quantified statistically. The 
phase mapping technique circumvents the difficulty of the circular characteristic of 0^ 
and illustrates the statistical significance of phase difference at the same time. This 
generalized method on phases allows us to detect weak coupling of phases on any Afc 
scales. 

Key words: methods: data analysis - techniques: image processing - large-scale 
structure of the Universe 



1 INTRODUCTION 

The morphology of the large-scale structure in the Universe 
is that of a complex hierarchy of nodes, filaments and sheets 
interlocking large voids. The Fourier-space description of 
such a pattern is dominated by the properties of the phases 
rather than the amplitudes of the Fourier modes (Chiang 
2001). According to the prevailing theoretical ideas this pat- 
tern developed by a process of gravitational instability from 
an amorphous pattern of density fluctuations characterized 
by a Gaussian field with random phases. Since the non- 
random phases of the present structure have grown from 
random-phase initial perturbations then there is strong mo- 
tivation for understanding how phase information develops 
within this paradigm and to construct a statistical descrip- 
tion of galaxy clustering that could be used as a test of the 
basic idea. 

Unfortunately, quantifying the properties of Fourier 
phases is difficult for a number of technical reasons, so 
their use in statistical studies has so far been limited. Ry- 
den & Gramann (1991), Soda & Suto (1992) and Jain & 
Bertschinger (1996) focused on the evolution of individual 
phases away from their initial values but since the initial 
phases are unknown these studies can not be used as the 
basis of a statistical descriptor. The pattern of association 
between phases is subtle and hard to visualize which makes 
a statistical test hard to construct a priori. 

As the first step in a different approach towards quan- 
tifying phase information, Coles & Chiang (2000) proposed 



a colour representation method to visualize phase coupling 
that at least reveals qualitatively how phase information 
arises during the evolution of A-body experiments but does 
not in itself constitute a statistical descriptor. In a related 
study, Chiang & Coles (2000) quantified phase information 
using a statistic derived from the Shannon entropy of the 
distribution of successive phase differences. This study dis- 
played interesting relationships between phase entropy and 
gravitational clustering but still did not provide a general 
statistical description. 

In this paper we use a generalization of the concept of 
a return map (May 1976; Chiang & Coles 2000) to transfer 
the phases of different Fourier modes on to a bounded square 
upon which simple statistical tests can be applied. In this 
way, we build upon the earlier studies (Chiang & Coles 2000; 
Coles & Chiang 2000) to construct a method that allows us 
to transform the phase information in a clustering pattern 
into a more useful form. 



2 PHASE COUPLING IN THE NONLINEAR 
REGIME 

The mathematical description of an inhomogeneous Uni- 
verse revolves around the dimensionless density contrast, 
8(x), which is obtained from the spatially- varying matter 
density p(x) via 
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8(x) 



p{x) - pQ 
Pa ' 



(1) 



where po is the global mean density. When the density per- 
turbation is small, the evolution of the density contrast can 
be obtained analytically through linear perturbation theory 
from 3 coupled partial differential equations. They are the 
linearized continuity equation, 

05 
dt 

the linearized Euler equation 



a 



(2) 



dv a 
at a 



— V x p- -V a 
pa a 



and the linearized Poisson equation 
V:t 2 ^ = 47rGa 2 po$. 



(3) 



(4) 



In these equations, a is the expansion factor, p is the 
pressure, Vz denotes a derivative with respect to the co- 
moving coordinates x, v — ax is the peculiar velocity and 
4>(x,t) is the peculiar gravitational potential. From Eq. (2)- 
(4), and if one ignores pressure forces, it is easy to obtain 
an equation for the evolution of 8: 

5 + 2{?-)6-4irGpo8 = 0. (5) 

For a spatially flat universe dominated by pressureless mat- 
ter, po(t) — l/6irGt 2 and Eq. (5) admits two linearly in- 
dependent power law solutions 5(x,t) — b±(t)8o(x), where 



8o(x) is the initial condition, b+(t) oc a(t) k t ' is the 
growing mode and b-(t) oc t^ 1 is the decaying mode. 

It is useful to expand the density contrast in Fourier 
series, in which 8 is treated as a superposition of plane waves: 



8(x) = 8(k) exp(ik ■ x) 



(6) 



The Fourier transform 8(k) is complex and therefore pos- 
sesses both amplitude |<5(fc)| and phase (j)^ where 



5(fe) = |5(fe)|exp(# fc ). 



(7) 



In the standard picture of 'gravitational instability' 
model for the origin of cosmic structure, particularly those 
involving inflation, the initial perturbations are Gaussian 
(Bardeen et al. 1986). The most relevant property of Gaus- 
sian random fields is that they possess Fourier modes whose 
real and imaginary parts are independently distributed. In 
other words, they have phase angles <\>y. that are indepen- 
dently distributed and uniformly random on the interval 
[0, 2n]. Terms in the perturbative evolution equations for 
the Fourier modes that represent coupling between different 
waves are of second (or higher) order in 8 and these are ne- 
glected in linear perturbation theory. When fluctuations are 
small, i.e., during the linear regime, the Fourier modes evolve 
independently (Eq.( 5) and (6)) and the Gaussian character 
is retained in the linear regime, where the phases remain 
independent and uniformly random. In the later stages of 
evolution, however, modes begin to couple together. In this 
non-linear regime that Fourier phases become non-random. 
For a thorough review of the theory and implications of non- 
linear evolution from the point of view of perturbation the- 
ory, see Bernardeau et al. (2002). 

Standard methods of analysis proceed via the power- 
spectrum, P(k), essentially proportional to |<5(fc)| 2 . The 




Figure 1. 2D TV-body simulations for initial power spectral in- 
dex, from left to right, n = —1, 0, 1 and 2, respectively. We 
choose 7 stages, named stage 1 to 7 from top to bottom, with 
an increasing level of non-linearity as described in the text. We 
control the simulations by setting the same initial random phase 
configuration for each index. Consequently, it is easy to see that 
the evolving structures have density concentration at the same 
locations; the difference is due to the initial spectral index: large 
n will produce more clumps, small n will have filaments. 



probabilistic properties of Gaussian random fields are com- 
pletely specified by knowledge of P{k). Higher-order quan- 
tities based on 8(k) can also be defined, such as the bis- 
pectrum (Peebles 1980; Matarrese, Verde & Heavens 1997; 
Scoccimarro et al. 1998; Scoccimarro, Couchman & Frie- 
man 1999; Verde et al. 2000; Verde et al. 2001; Verde et 
al. 2002), which vanishes for Gaussian fields, or quantities 
related to correlations of |<5(fc)| 2 (Stirling & Peacock 1996). 
Phase coupling results in a non-Gaussian field in which the 
bispectrum and higher-order polyspectra may be non-zero 
(Watts & Coles 2002). Phase information is at the heart of 
non-linear galaxy clustering. 



3 DIRECTIONAL PHASE MAPPING 

There are two principal difficulties involved in constructing 
a statistic from Fourier phases. One is that because phases 



© 0000 RAS, MNRAS 000, 000-000 



Return mapping of phases and the analysis of the gravitational clustering hierarchy 3 



stage 3, n='1 at (rn i n)=(O s 33) 




Random phc 



it (m,n) = (-2,44) 



1 00 200 300 




(a) 



(b) 



stage 2, n=1 at (m,n) = (-2,44) 



stage 3, n=1 at (m,n) = ( 



300 



200 



100 




300 



200 




(c) 



(d) 



Figure 2. Phase mapping onto the return map. The horizontal and vertical axes represent (pj, and </>fc +/ \fc, ranging from 0° to 360°. In 
panel (a) the phases arc taken from the stage 3, spectral index n = 1 of the .W-body simulations shown in Fig. 1 and the scale of phase 
mapping is Afc = (m, n) = (0, 33). The phases are condensed along the diagonal strip. Panels (b), (c) and (d) are the smoothed versions 
of various phase maps where the contour levels are drawn upwards starting from the mean value. Panel (b) is from a random-phase 
realization and (c) is that of stage 2, n = 1. For comparison, both (b) and (c) are mapping from the same Afc scale (m,n) = (—2,44). 
Panel (d) is the smoothed version of panel (a). 



reflect the morphology, their values change according to the 
position of structural features (Chiang 2001). For exam- 
ple, the phases of a peak in the form of Dirac-5 function 
Sd(x — xo), 4>k = kxo suffer change in slope along the fc-axis 
when there is shift of the peak 5d(x — xq — x ), the phases be- 
ing (f>k = k(xt)-\-x ). If a pattern is statistically homogeneous, 
any descriptor of it should be translation-invariant and this 



is manifestly not true of the phases themselves. The other 
problem is that phases are of circular measures and there- 
fore defined modulo 2tt. Traditional measures of association, 
such as covariances of the form (4>i4>j), are based on the as- 
sumption that the measure associated with the variable is 
linear and are therefore not appropriate to cases where val- 
ues separated by 2ir are in fact equal in value. 
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To address these problems, Chiang & Coles (2000) used 
the phase difference (or phase gradient) , Dk , defined in one 
dimension by 

Dk = 4>k+\ — 4>k, (8) 

i.e. for neighbouring phases. In two or three dimensions dif- 
ferences can be taken in orthogonal directions. The quantity 
Dk has the twin advantages that for random phases it is also 
random but upon translation x' it changes by a constant x' 
for all k. The statistical properties of the set of differences Dk 
contain information about the correlations of neighbouring 
Fourier modes. Strong correlation of the neighbouring modes 
at large k corresponds to the highest peak in the clustering 
pattern. Naselsky, Novikov and Silk (2002) used this charac- 
teristic to extract point sources in the CMB map. Moreover, 
Chiang et al. (2001) used the phase analysis for extracting 
the in-flight beam shape properties of CMB experiments. To 
construct a more general description of phase coupling we 
need to extend this method to modes that are not necessarily 
neighbours. We do this by constructing a directional phase 
map, based on the return maps used in non-linear dynamics 
(May 1976). 

The basic idea is simple and based on a study of one- 
dimensional examples contained in Chiang & Coles (2000) 
which provides a useful illustration of the more general ap- 
proach. With a set of phases (f>k from the Fourier trans- 
form of a one-dimensional process, one can plot a map of 
4>k against <pk+i f° r each pair (<f>k,4>k+i)- If the phases are 
random this will be a scatter plot with points distributed 
randomly within the bounded square of side [0, 2ii\. If there 
is association between neighbouring phases the plot will con- 
tain correlations; the quantity Dk is sensitive to linear as- 
sociation. If the spatial pattern consists of a single high- 
amplitude peak the points display linear association and are 
mapped into a diagonal lines on the diagram. 

In what follows, for illustration, we shall use two- 
dimensional examples based on numerical simulations with 
periodic boundaries, so we take Afc = (m, n) where m and 
n are integers. The simulations are done on a 512 2 grid. In 
Fig. 1 we show 4 sets of such iV-body simulations for initial 
power spectral index n = — 1, 0, 1 and 2 (see Chiang & Coles 
2000 for details of the simulations). The evolutionary stages 
are characterized by an increasing scale of non-linearity de- 
fined by ((8p/p) 2 ) k = b\(t) J () feNL P(k) d 2 k = 1, where 

b + (t) is the growing mode of the linear density contrast and 
P(k) is the linear extrapolation of the initial power spec- 
trum. This definition of &nl identifies the corresponding 
scale 27rfe^ L as the boundary between linearity and non- 
linearity. The stages in Fig. 1 are chosen such that the scales 
of non-linearity ^nl between any two successive stages vary 
by a factor of 2. The levels of non-linearity of the stages are 
thus fc NL = 256, 128, 64, 32, 16, 8 and 4* 

We can use these simulations to illustrate how we 
extend the return mapping between neighbouring phases 
(Afc = 1) to pairs of phases with any Afc scales in fc-space. 
We map all pairs of phases 4>(i,j) and <j>(i + m,j + n) onto 
the x and y values of the return map. The axes therefore 

* To avoid confusion with the panels in the captions, we re-name 
the stages as 1-7, which are originally named as stage a-g in Chi- 
ang & Coles (2000). 



range over [0, 27r] for both 0(fe) (x) and 0(fc + Afc) (y) axes. 
For example, for (m, n) = (4, 6) we have points on the re- 
turn map (<t>(i,j),<f>(i + 4, j + 6)) for all i G [-255,256], 
j € [1, 256], i.e., all points (0(1, 1), 0(5, 7)), (0(2, 1), 0(6, 7)), 
(0(1, 2), 0(5, 8)) . . ., from a 2D Fourier transform of a real- 
isation. This represents the directional phase coupling for 
coupling scale Afc = (Ak x , Ak y ) = (m,n) in fc-space. The 
neighbouring phase differences in the fc^-direction and k y - 
direction used by Chiang & Coles (2000) simply corresponds 
to (m,n) = (1,0) and (0, 1) respectively. 

In Fig. 2 (a) we show one example of phase mapping 
from the realization of stage 3, spectral index n = 1 simula- 
tion. The particular Afc in this example is (m,n) = (0,33). 
This panel demonstrates how weak coupling between phase 
pairs with fixed scale Afc can manifest itself in the return 
map as non-uniform density in the map plane. 

This directional phase mapping approach circumvents 
the problem of the circular character of 0^ but does not at- 
tempt to condense all the related information into a single 
quantity. It, on the other hand, exploits all the information 
between all Fourier modes. For example, the phase coupling 
of a ID distribution can be expressed in a (2D) return map. 
It therefore allows us to build simple statistics to test the 
significance of general non-randomness. The phase difference 
between any pair at a fixed scale becomes a single point on 
the return map for that scale. The circular characteristic of 
0^, is transferred to a bounded square, topologically equiv- 
alent to a torus owing to the periodicity of x and y axes. 
The bands seen in Fig. 2 (a) therefore correspond to twisted 
linear features on this torus. 

The key advantage of directional phase mapping is that, 
for a Gaussian random field, any directional phase mapping 
for any scale Afc = (m, n) should produce a random Poisson 
distribution. Weak phase coupling will produce correlations 
at large vectors (m, n) while strong non-linearity will pro- 
duce highly non-uniform phase maps at all scales (m,n). 



4 A x 2 TEST ON PHASE MAPS 

Once we have transferred the phase information onto a phase 
map like that shown in Fig. 2 (a), many different statistical 
tests can be used to analyse its properties. Here we outline 
a simple yet powerful method. 

First we smooth the return map. Smoothing enhances 
our visualization of the pattern of phase coupling. In Fig. 2 
(b), (c), and (d) we show the contour maps of the smoothed 
return map. In these simple illustrative experiments we di- 
vide the square of the return map into 128 2 pixels, and we 
bin the 32 768 points so that, for a perfectly even distribution 
in the map plane the occupation of each pixel is {p(i, j)) = 2. 
Then we smooth this 128 2 mesh by 




The contour levels are drawn upwards starting from the 
mean value. Panel (b) is the contour of a smoothed return 
map from a realization of random phases. Panel (c) is that 
from stage 2, n = 1 of Fig. 1 with (m,n) = (—2,44). For 
comparison, the coupling scale Afc for both (b) and (c) is 
set the same. We can see that even for the mildly non-linear 
regime represented by stage 2 (for which 5 ~ 1.8), phase 



© 0000 RAS, MNRAS 000, 000-000 



Return mapping of phases and the analysis of the gravitational clustering hierarchy 5 




Figure 3. The x 2 statistics on grey scale for different realisations on the (m, n) plane. The 4 panels correspond to the 4 realisations of 
Fig. 1, stage 2, 3, 4 and 5, spectral index n = 1 of the 2D iV-body simulations. In order to show clearly the maximal points, we set the 
same grey scale for panel (a), (b) and (c), and note that the grey scale of (d) is different. The maximal x 2 for stage 2 is 2.53 X 10 — 2 at 
(m,n) = (-2,44), and 3.27 X 10~ 2 for stage 3 at (m,n) = (0,33), both mapping of the particular scale being shown in Fig. 2 (c) and 
(d). For panel (c) and (d) the X 2 m ax = 4 -90 x i0 ~ 2 and !- i7 x at (.m,n) = (12,5) and (-4,8), respectively, which indicates the 

scale of phase coupling in the highly non-linear regime is at small Afc. 



mapping does reveal the existence of coupling by starting to 
condense on the diagonal strip. Panel (d) is the smoothed 
version of (a) with (m,n) = (0, 33), in which the condensed 
strip from (a) is much clearer. 

We define a mean x 2 statistic as 



?= ^ E kM^ (10) 

where M is the number of pixels we assign on the return 
map and p is the mean value for each pixel. 

In Fig. 2, the smoothing scale on the 128 2 mesh is R — 2 
and the contour levels are drawn starting from the mean 
value. We use this pixel size and smoothing scale in all our 
calculations of x 2 in the following, but one can vary the scale 
as part of a statistical test. 



5 SIMULATION RESULTS 

We now illustrate the results of this analysis using the 2D N- 
body simulations described earlier. First we Fourier trans- 
form the realizations of the N 2 mesh (N = 512 in our simula- 
tions). Because of the reality of the original distribution, and 
the consequent Hermitian conjugate relations in the Fourier 
image, only half of the Fourier transform contains indepen- 
dent information. We end up with N 2 /2 Fourier modes avail- 
able, so we take —k N / 2 + 1 to k N / 2 in k x axis and from 1 to 
fcjv/2 in k y axis. 

We carry out a phase mapping for each (m, n). The di- 
rectional phase mapping is performed for the vectors (to, n), 
where to £ [—49,50] and n £ [1,50] for the phase base 
4>(i,j), where i £ [-127,128], and j £ [1,128], a quarter 
of the available phases. The limited range of mapping vec- 
tors is chosen to ensure that the map can be constructed 
without running out of sensible wavenumbers. Thus, in such 
a case, for each return map of (m, n) there are 2 x 128 2 
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Figure 4. The maximal \ 2 ls plotted against ID Fourier scale 
|Afc| = Vm 2 + n' 2 for stage 2 and 3 of simulations from initial 
spectral index n = — 1, 0, 1 and 2. The maxima are taken from 
each |Afc| ring of the (m, n) plane. 



Figure 5. The maximal x 2 ls plotted against ID Fourier scales 
|Afc| = V rn' 2 + n' 2 for stage 4 and 5 of simulations from initial 
spectral index n = — 1, 0, 1 and 2. The maxima are taken from 
each |Afc| ring of the (m,n) plane. 



points. We set up M = 128 2 pixels in this [0, 2-k] square and 
then smooth the return map to decrease the fluctuation. The 
X 2 (w-, n) for each fixed scale (m, n) is calculated by Eq. (10). 

Figure 3 shows 'supermaps' of all the phase information 
contained in the x 2 { m ,n) using a grey scale in the (m,n) 
plane from various realisations of simulations. We are inter- 
ested specifically on the mild non-linear regime, so in Fig. 3 
the panels (a), (b), (c) and (d) correspond to stages 2, 3, 
4 and 5 of n = 1 in Fig. 1, respectively. The m-axis ranges 
from —49 to 50 and n-axis ranges from 1 to 50. Each pixel 
therefore displays the level of phase coupling on a certain 
scale (Ak x , Ak y ) = (m, n) in terms of the \ 2 statistics. The 
bright points in Figure 3 are direct indications of phase cou- 
pling for the corresponding Afc scales. 

The reason we present the specific scales (m, n) = 
(—2, 44) for stage 2 and (m,n) = (0, 33) for stage 3, n = 1 in 
Fig. 2 is that the calculation of the \ 2 statistics shows them 
to be 'hotspots' of phase correlation. The scales of phase 
coupling for those two panels correspond to the brightest 
points, i.e. the maximum \ 2 on the corresponding (m,n) 
planes in Fig. 3. 

In Fig. 4, 5 and 6 we plot the maximum \ 2 against Afc 
for all 4 sets of simulations. These ID plots show that the 
maximal phase coupling from each ring | Afc| = V m 2 + n 2 of 
the (m,n) plane. As we control the simulations by assigning 
the same set of initial random phases, These ID plots will 
display how the scale |Afc| of phase coupling is related to 
morphology. 



The iV-body simulations we have carried out is of self- 
similar nature, that is, a distribution function f(x,t) has 
the same statistical measure as the re-scaled one 

/ i-> ^fix/X 13 , \t), as t h-> At. (11) 

With the reciprocal-scaling property of Fourier transform, 

'<~> = R F (I)' (12) 

for non-linear scales xnl increasing, i.e. a > 1, the corre- 
sponding scales in k space are decreasing. Although phases 
do not possess a linear relationship owing to their circular 
nature, it can be understood qualitatively that, if there ex- 
ists coupling between pairs of phases with fixed Afc, this 
scale has to decrease as gravitational clustering proceeds. 

It is therefore clear that for the case such as n = — 1, 
where large-scale filaments are the prominent feature, phase 
coupling starts from low Afc, which also has cascade effect 
on to higher Afc, as phases strongly couple on any Afc might 
also do so at multiples of Afc. For high n, on the other hand, 
where small clumps form first, phase coupling starts from 
large Afc. This also explains why the Shannon entropy from 
neighbouring phase difference can produce considerable re- 
sults at early stages for n = — 1, but not for n = 2 (see the 
fig. 5 of Chiang & Coles 2000). On the other hand, the cou- 
pling between the amplitudes are enhanced by the factor 1 /a 
as clustering proceeds, so mode-mode coupling between am- 
plitudes at early stages is not as obvious as between phases. 
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Figure 6. The maximal x 2 ls plotted against ID Fourier scales 
|Afc| = V rn' 2 + n' 2 for stage 6 and 7 of simulations from initial 
spectral index n = — 1, 0, 1 and 2. The maxima are taken from 
each |Afc| ring of the (m,n) plane. Note that the scale of j/-axis 
is 10 times larger than that in Fig. 4 and 5. 



For a discussion of how this relates to the development of 
phase correlations in mildly non-linear evolution, see Watts 
& Coles (2002). 

The x 2 statistics shown on the (m, n) plane and those 
ID plots confirm the visualization of phase coupling pre- 
sented by Coles & Chiang (2000). Phase coupling firstly ap- 
pears on large Afc when the the scale of non-linearity is small 
in real space, then it shifts on the (m, n) plane to small Afc 
in Fig. 3 (c) and (d) and finally dominates at neighbouring 
modes as seen in Fig. 6. 



6 CONCLUSION 

We have generalized a method based on phase mapping on 
the return map. This simple, easy-to-implement method can 
detect phase coupling at any scales Afc in k space. We apply 
this method to two-dimensional simulations of gravitational 
clustering and the result has shown that even when the evo- 
lution is in the mild non-linear regime, phase coupling on 
certain scale is revealed through the \ 2 statistics on the 
(m, n) plane. 

In contrast to other methods, such as the Shannon en- 
tropy of the distribution of neighbouring phase differences 
(Chiang & Coles 2000), this method does not require large 
number of phases. Moreover, this approach can detect the 



scale of phase coupling through the phase mapping as shown 
in Fig. 2. 

With the systematic iV-body simulations shown in 
Fig. 1, we have also demonstrated in Fig. 4-6 that the scale 
of phase coupling differs according to the clustering mor- 
phology: modes between small Afc for large-scale filaments, 
large for small clumps. 

This method reveals a signature of non-linear gravita- 
tional instability, but also offers the opportunity to provide a 
general test of Gaussianity that could be applied to cosmic 
microwave background temperature maps. In future work 
we shall evaluate the effectiveness for such method. 
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